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Abstract 

The Squared- Gradient approximation to the Modified-Core Van der Waals density functional 
theory model is developed. A simple, explicit expression for the SGA coefficient involving only 
the bulk equation of state and the interaction potential is given. The model is solved for planar 
interfaces and spherical clusters and is shown to be quantitatively accurate in comparisons to 
computer simulations. An approximate technique for solving the SGA based on piecewise-linear 
density profiles is introduced and is shown to give reasonable zeroth-order approximations to the 
numerical solution of the model. The piecewise-linear models of spherical clusters are shown to be 
a natural extension of Classical Nucleation Theory and serve to clarify some of the non-classical 
effects previously observed in liquid-vapor nucleation. Nucleation pathways are investigated using 
both constrained energy-minimization and steepest-descent techniques. 
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I. INTRODUCTION 



The description of inhomogeneous systems remains one of the most important problems 
in several areas of physics. Recently, it has been shown that classical Density Functional 
Theory (DFT) can give quantitatively accurate results for many inhomogeneous systems 
including the structure and surface tension of the liquid-vapor interface, confined fluids near 
walls, in slit pores[l| and in small cavitiesj2| and of liquid- vapor nucleation sl- While these 
results confirm that DFT can be a useful tool for making quantitatively accurate calculations, 
they provide limited physical insight due to the complexity of the DFT models and the need 
for extensive numerical calculations. 

One simplification of DFT is to relate it to more physically-inspired descriptions of in- 
homogeneous systems. The oldest, and in many contexts most important, such model is 



the squared-gradient approximation (SGA) that dates back to van der Waals[ 



HE 



5| and has 



been re-invented and exploited in many different circumstances, most notably by Landau in 
the context of phase transitions!^ and by Cahn and Hilliard[7| for the description of inter- 
faces. Evans established the link whereby SGA is an approximation to more fundamental 
DFT Is I and this link has more recently been refined and systematized 9-11 1. Despite this 
justification for SGA, it nevertheless has the reputation of giving only a qualitative, at best 
semi-quantitative, approximation to real systems 12|. One result reported here is that SGA 
based on DFT can in fact be surprisingly accurate. 

If DFT can be well approximated by SGA, then a substantial reduction of complexity 
has been achieved. However, for most applications, SGA must also be solved numerically. 
No matter what details go into the models, this fundamentally involves solving for the 
spatially-varying density which minimizes the free energy. A second development presented 
here is the use of piecewise-linear approximations for the density profiles. In the simplest 
case, this involves approximating the interface between two bulk phases as a linear function. 
The calculation can be systematically improved by replacing the single linear function by 
a sequence of linear functions, or links, so that the true profile is obtained in the limit 
of the number of links going to infinite. It is shown that the simplest single-link profile 
provides a rather accurate approximation to the numerical solution of the SGA with the 
benefit of giving nearly analytic expressions for surface tensions and interfacial width of 
planar interfaces and of interfacial width, surface tension and Tolman length for spherical 



clusters. An additional benefit is that by reducing the DFT integral theory to a simple 
algebraic model, it gives a systematic link between DFT and Classical Nucleation Theory 
(CNT). 

Finally, the algebraic model for spherical clusters is used to study the nucleation pathway 
for homogeneous liquid-vapor nucleation. A first method is to minimize the free energy 
subject to constraints such as fixed number of atoms in the cluster or fixed radius of the 
cluster. It is found these and other constraints are all problematic and fail to give a consistent 
picture of homogeneous nucleation. An alternative method is the construction of steepest- 
descent paths in density space linking the transition state (i.e. the critical cluster) to the 
bulk liquid and bulk vapor free energy minima 13|]- This is found to give a simple and 
plausible description of the nucleation pathway. 

In the next Section, the SGA is formulated based on the quantitatively accurate Modified- 
Core Van der Waals model DFT. The result is an SGA model for fluids that requires only 
the bulk equation of state and the interaction potential as input. To the extent that the 
equation of state can be accurately approximation by, say, thermodynamic perturbation 
theory, only the interaction potential need be specified. The application to planar interfaces 
and for spherical clusters and the piecewise-linear approximation are also discussed. In the 
third Section, this model is used to explore liquid-vapor nucleation. First, the connection 
with the ideas that underlie Classical Nucleation Theory is made. There, different methods 
of formulating the problem of the determination of the nucleation pathway are described and 
their relative merits compared. The fourth Section gives a comparison of the SGA and the 
analytic models to computer simulation of both planar interfaces and clusters. The paper 
ends with a discussion of the results. 



II. THEORY 



A. Squared-gradient approximation 



Density Functional Theory is an approach to equilibrium statistical mechanics which is 



formulated in the grand-canonical ensemble at constant temperature, T, chemical 



otential. 



B Q,li5i. 



It 



/i, and volume V. The fundamental object is the local number density p (r) 
can be shown that there exists a functional of the local density, Q [p], having the property 



that it is minimized by the equihbrium density function and that its value at this minimum 
is the grand-potential for the system. It can be written as 

^[p] = F[p] + J (0 (r) - /i) p (r) rfr (1) 

where (j) (r) is any external field that may act on the system and where F [p] is independent 
of the field but otherwise unknown except for special cases such as the ideal gas for which 

Fid [p] = kBT j p (r) (Inp (r) - 1) dv (2) 

where is Boltzmann's constant. When evaluated at the equilibrium density, F [p] is the 
intrinsic Helmholtz free energy. Minimization of Vt [p] gives the Euler-Lagrange equation 

f(^ + ^W-. = (3) 

A key result used to guide the formation of models is that the excess functional Fex [p] = 
F [p] — Fid [p] is related to the two-body direct correlation function via 

6^I3F^^ [p] 

c ri, rs; p = . ...... (4) 

6p{ri)6p{r2) 

where /3 = l/ksT. By definition, a bulk fluid is a system with constant equilibrium density, 
p (r) = p in which case F [p] = F (p) is the usual Helmholtz free energy function for the 
fluid. 

It can be shown llj that a systematic expansion of F [p] takes the form 

F[p]= |^{/(p(r)) + iA^(p(r))(Vp(r)f + ...}rfr (5) 

where / (p) = yF (p) is the Helmholtz free energy per unit volume of the homogeneous 
system, the ellipses indicate higher order terms in the gradients and the coefficient of the 
second order term is 

{3K (p(r)) = -^jc {r^2■, p(r)) rl,dr^dr2 (6) 

where c(ri2;p) = c(ri,r2;p) is the translationally invariant direct correlation function of 
the bulk state. The gradient expansion is thus seen as a link between the properties of the 
bulk state - which are generally accessible - and those of the inhomogeneous state which 
are generally much more difficult to determine. Truncating the expansion at second order 



gives the Squared- Gradient Approximation (SGA). Substituting into Eq.([3]) and taking the 
external potential to be zero gives 

V-Kip (r)) Vp (r) - i (^A' (p (r))) (Vp (r))^ - ^^^M^ = (7) 

where u (p) = f (p) — pp is the grand potential per unit volume. 

To implement the theory, two elements are necessary: the Helmholtz free energy and the 
direct correlation function of the bulk system. These are not independent: the equation 
of state is easily calculated from the DCF via the compressibility equation as discussed 
below. For realistic potentials, the DCF can be calculated using liquid state theory such as 
the Percus-Yevick or the HNC approximations. However, when dealing with an interfacial 
system, the density typically varies a lot: for a liquid- vapor interface it obviously spans the 
range from (dense) liquid to (low density) vapor. In particular, it passes through densities 
that lie in the two-phase region of the bulk phase diagram and liquid-state theories often 
do not have solutions in those regions. While methods of circumventing this problem by 
means of e.g. interpolation have been proposed (see e.g. Ref. 12|) this remains a problematic 
issue. Of course, this is an issue facing all DFT's since the free energy functional is always 
related to the DCF so that it might be suspected that a successful DCF would imply some 
means around this problem. The Modified-Core Van der Waals model DFT is based on a 
simple approximation to the DCF and gives good results for a wide variety of interfacial 
systems [ll,|3|. The idea behind it is to begin with the simplest Van der Waals model whereby 
the DCF is approximated as that of a hard-sphere system with a mean-field treatment of 
the attractive tail of the interaction, 

cvDW (ri2; p) = chs (ri2; P, d) - (3v (ria) Q (r - d) (8) 

where v{r) is the molecular pair interaction potential, d is the effective hard-sphere diameter, 
Chs (ri, p, d) is the hard-sphere DCF and 0(x) = 1 for x > and zero otherwise. Because 
of the relation between the DCF and the excess free energy of the uniform bulk fluid, 

Pfexip) = Pf ip) - f3fid{p) = - y J dp2 J J J c(ri2;pi)rfridr2, (9) 

this implies a rather inaccurate equation of state. Furthermore, it is discontinuous at the 
hard-sphere boundary. To address both of these problems, a linear correction is added to 
the core region giving the full approximation 16 1 

c (ri2; p) = CHS (ri2; P, d) + (bo + 6 (rf - r) - (3v (ru) Q {r - d) (10) 



The coefficients bo and bi depend on both density and temperature and are chosen to re- 
produce a given equation of state and to make the DCF continuous. This has been shown 



to give a rather good approximation for dense fluids 16|. At low density, one has the exact 
result 

lim c (ri2; p) = - (l - e-^^'^^^^^) (11) 

and is clearly not reproduced by the simple model. This limit can be enforced by generalizing 
to 

c (ri2; p) = CHS (ri2; P, d) + (c (ria; 0) - chs [ru] 0)) + (bo + 6 (rf - r) (12) 

which also gives a reasonable approximation. However, as most properties are insensitive to 
the low-density limit of the DCF, it is more convenient and not much less accurate to work 
with the simpler approximation given above. 

Assuming the hard-core diameter is chosen to be independent of density, the core- 
correction coefficients are fixed by the requirements that the model reproduce the given 
bulk equation of state, / (p), and that the DCF be continuous giving 

(p) = ^fHS (p) - ^TTd' (-bo + -b^\ + Arc fiv (r) r^dr (13) 
Chs {d-] p, d) + bo + bi = -(3v (d) 

As shown in Appendix [Al using common models for the hard-sphere DCF, the final result 
for the SGA coefficient takes the simple form 

nd^ ... d^ , , 2n , . „ dBv (r) , , 

where {'pd^) depends on the particular hard-sphere model used (explicit expressions are 
given in Appendix |A]) . However, since the coefficient of a^ip) is so small, this term con- 
tributes only about 1% to the over-all value of K and can usually be safely neglected. Thus, 
in this model, the SGA coefficient is a relatively simple function of the hard-sphere diame- 
ter, the potential and the equation of state. While any reasonable value for the hard-sphere 
diameter could be used, the calculations presented below are based on the Barker-Henderson 
formula Q, El, 

d= {l-exp{-/3v{r)))dr, (15) 
Jo 

where Tq is the smallest solution of v (ro) = 0. 



B. Planar and spherical interfaces 



A stable planar interface can only exist when the liquid and vapor are at conditions 
of coexistence so that the value of the chemical potential is ficoex and the bulk liquid and 
vapor densities are pi^coex and Pv,coex respectively. (All of these quantities are temperature- 
dependent.) Then, one can impose a density profile that varies in only one dimension, p{z) 
so that Eq.Q becomes 

'-K (p (.» A, (,) _ 1 f (p (.))) ( ±p (.)) " - '-^MeH .0 (16) 



dz dz 2 \dp{z) J \dz J dp{z) 

Assuming that the density takes the value of the bulk liquid and vapor far from the boundary, 
and noting that u {pi,coex) = ^ {Pc,coex) = ^coex the profile equation can be integrated to get 

K{p{z))^—p{z)j =2(a;(p(^))-a;,,,.) (17) 

This allows the excess free energy, hereafter referred to as the "surface tension", to be 
evaluated as 

n-n. 



A 



= r'"°'' V2 {oo (p) - uj,oex) K {p)dp (18) 

while the actual profile must be obtained by integrating Eq. (lT7|) numerically. 

Away from coexistence, a planar interface is unstable. In this case, it is more pertinent 
to study clusters since an unstable phase will transform to a stable phase by the formation 
and subsequent growth of a critical cluster. Assuming spherical symmetry, the equation for 
a (met a-) stable profile becomes 

K (P (r)) --rr.rp (r) + " (0 - ttTTIV^^ (p (^)) hlP (0 = 



r dr"^ dr dr 2\dp{r) J \dr J dp{r) 

(19) 

and this does not admit of an exact quadrature. The only nontrivial solution will correspond 
to the critical cluster in which the density near the origin will be close to that of the stable 
phase, while the bulk (i.e. the density far from the origin) will be that of the unstable phase. 

C. Analytic approximations 

An alternative to solving these equations numerically is to assume some ansatz for the 
density, p (r) = p (r; F), where the quantity on the right has a specified spatial dependence. 



e.g. a sigmoidal function for the case of a planar interface, and where F represents a collection 



of parameters (e.g. the width and center of the sigmoidal function). In 



act, this procedure 



was been used for the planar interface by Telo da Gama and Evans 1 181] and, recently, for 
droplets by Ghosh and GhoshpJ?]. Then, rather than solving Eq.([7]) to get the profile, 
one substitutes the ansatz into the expression for Q [p] and extremizes the free energy with 
respect to the parameters F. For most cases, this will still involve numerical calculations 
however more progress is possible if one assumes the simplest reasonable approximation 



which is a piecewise-continuous profile. Thus, for the planar profile, one can try 



P—ooi Z 2 

P-oo + (Poo - P-oo) -f < ^ < f (20) 



2 ^ ^ 2 

Poo) "2" ^ ^ 

There are three parameters characterizing this profile: the densities p±oo on either side of the 

interface and the width of the interface. The profile is continuous and differentiable except 

at z = ±Y but this is sufficient to allow evaluation of the free energy. It is easy to see that 

the minimization of the free energy for the case that the volume is very large requires that 

df (Poo) _ dfip^oo) _ 
dpoo op-00 

so that a non-trivial interface is only possible at coexistence in which case one of the densities 
must be that of the coexisting liquid and the other that of the coexisting vapor. Then, the 
excess free energy per unit area is 

7 = ^^^= r U (P (^)) - ^c.ex. + Ik {P (z)) ( ^°° "J^'^ y] (22) 



or, more simply, 



with 



^ = w{uJ- a;,_) + )V (23) 



2 /"Poo 

u = / u (p) dp (24) 

Poo ~ P-oo J 

Poo 



K = — / K (p) dp 

'p-oo 



1 

Poo P— 00 

(For simplicity of notation, the dependence of u and K on the densities is not indicated 
explicitly). Minimizing with respect to the width gives 



and 

7 = (poo - P~oof ip - Wcoex) K (26) 

which is very similar to the exact result given in Eq.( !T8|) . However, in the present case, 
one has analytic expressions as well for the entire profile, including the width which can be 

written as 

w = --^ ^ ^P--P-)"K (27) 

2[UJ - Ucoex) 7 

giving a simple relation between the interfacial width, the SGA coefficient and the surface 
tension. Notice that using the explicit expression for the SGA coefficient and neglecting the 
small hard-sphere term, one has that 

since the integral of the density-dependent part of the coefficient gives no contribution at 
coexistence. In fact, if the density dependence of K{p) is weak, as will be seen to be the case 
for a simple fluid, then the approximation K{p) — )■ K should be adequate thus giving an 
even simpler expression for the gradient coefficient based solely on the interaction potential. 
One can make a similar ansatz for the spherical cluster. 



p[r) 



po, r < R 

Po + iPoo-po)"^, R<r<R + w (29) 
Poo, R + w < r 



so that, with Q^o = Vu{poo) and AQ = Q — Qoo, 

An = —R^Aoo (30) 

.,.n^. ( . ... q) ... . (i-^) (7r„ . Q . 

where V is the total volume, the density moments of the bulk free energy and SGA coefficient 
are defined as 

UJn (Poo, Po) = / (W (P) - (Poo)) dp (31) 

Poo-Pojpo VPoo-Po/ 

^ P - Po 



Kn (Poo, Po) = / K (p) dp 

poo -pa J \Poo-poJ 

Note that the density arguments have been suppressed (i.e. cJo = Auq (pocPo), etc.), that 
Au = (po) — uj{poo), and that the zeroth-order moment uJo{poo,Po) is the same as the 

9 



quantity cJ(poo,Po) characterizing the planar interface. Equation (130!) is completely general 
and involves no assumptions regarding the size of the cluster. 

As shown in AppendixlEl these piecewise-linear approximations can be systematically 
extended to include an arbitrary number of linear "links" in the profile, each with a sep- 
arate width and slope. It is expected that as the number of links grows, the profile will 
become closer and closer to the true, free-energy minimizing profile so that this provides 
an alternative form of solution of the unconstrained problem. As will be noted below, only 
a few links are needed to obtain high accuracy in the excess free energy. This therefore 
provides a systematic, and computationally cheap, alternative to the direct solution of the 
Euler-Lagrange equation for the profile. 

III. LIQUID- VAPOR NUCLEATION: EXTENDING CLASSICAL NUCLEATION 
THEORY 

For a given temperature, there is a unique value of the chemical potential at which the 
liquid and vapor phases can coexist. Any other value of the chemical potential implies 
stability of one phase over that of the other. In this Section, the transformation from the 
metastable phase to the stable phase is discussed based on the piecewise-linear model for 
the spherical cluster. The reason for concentrating on this model, rather than solving the 
SGA equations numerically, is to make contact with Classical Nucleation Theory and to 
show how it can be extended in a natural way to include non-classical effects. Furthermore, 
this simplified version of the theory is expected to be useful as a starting point for studying 
more complex systems. 

A. Classical Nucleation Theory 

For a given value of /i 7^ ficoex there is a metastable phase, pm, and a stable phase, ps 
and u (pm) > <^ (Ps) while u' (p^) = (p^) = 0. The system is initially in the metastable 
state p (r) = p.^ and homogeneous nucleation proceeds by the spontaneous formation of a 
critical cluster. Throughout this discussion, it will be assumed that such clusters are always 
spherical. In CNT, the interface between the two phases is sharp so that the cluster has a 
well-defined radius, R and, on physical grounds, the free energy of a cluster is assumed to 

10 



be given by 

^CNT = yi2'cu (Po) +\V- —R'j tu (poc) + AttR^ (32) 

where 7c is the planar excess free energy per unit area at coexistence and where the internal 
and external densities, po and poo are to be determined. Extremizing the free energy gives 
(Po) = (Poo) = showing that the interior and exterior densities are those of the bulk 
stable and metastable phases respectively. Then, since u (p^) < u (pm) and 7c > 0, the free 
energy necessarily has a maximum at the critical radius, 

UJ (pm) - W (ps) 

giving the free energy barrier 

IGvr 7"^ 

^CNT - ^ (pm) = , " , ,,2 - (34) 

[UJ {pm) - [ps)) 



This can then be used to estimate the nucleation rate[20| where, in the course of the analysis, 
Eq.( l32|) is used to determine the number of non-equilibrium clusters of a given size under the 
assumption that the number of clusters of size = ^R^po is proportional to the Boltzmann 
factor, exp {—/SQcnt (N)). Since there is only one variable that can vary between clusters, 
i.e. the radius R, the implied nucleation pathway is simply one of increasing radius. 



B. Energy- minimized pathways 



Comparison of Eq. (l32l) and (l30l) shows that the CNT expression results from the DFT 



model if one takes the limit ^ with 



Ko held fixed which suggests looking at 



the large R behavior of the DFT model. Suppose that the width in the planar (i.e. large 
radius) limit is Wq. Then, one can obtain the large-cluster limit by expanding in the small 
parameter e = Wq/R so that 



w = wq + ewi + ... 
Po = Poo + epoi + ••• 



(35) 



Minimizing the excess free energy with respect to w and po at fixed R then gives 



An = y {u (poo) - UJ (poo)) + AnR^-f 



-Wq 



Wo 



(36) 
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with the zeroth and first order densities determined from 

duj (poo) 



(37) 



00 



3 Wo / Aw (poo) , K (poo) 

POI = —-777-^ + 



dp. 

Poo - Poo UJ (poo) V ^0 Kq 

where it is understood that ujq = ujo {Poo, Poo): etc. The first equation specifies that the 
density in the large R hmit is that of the bulk liquid for the applied chemical potential as 
expected. The zeroth order width is found to be 



/ (POO — Poaf^ /na\ 

and the coefficients for the expansion of the surface tension term are 

7 = 2tf;oAa;o (39) 
uJi Ki 1 poi [u (poo) - ^ {poo) , K (poo) 

= h =^ + h 



Wo Kq 4 Poo - Poo V ^0 Kq 

Although the radius parameter i? is a model-dependent quantity, it can be related to a more 
physical quantity, namely the equimolar radius Re, via 

Rl = / (P (r) - p (oo)) dr = R' + \w {6R^ + ARw + w^) . (40) 

While this simple calculation serves to establish a direct link between the DFT and the 
CNT free energy model, it leaves open the question of what the nucleation pathway might be. 
There is certainly no reason to assume that clusters grow by increasing the radius parameter 
in this model while minimizing the free energy with respect to the other parameters. In fact, 
as shown below, there is reason to believe that the radius parameter varies non-monotonically 
along the nucleation pathway and that it cannot be used to parameterize movement along 
the nucleation pathway: i.e. it is not a good reaction coordinate. Physically, one expects the 
cluster to grow by adding atoms so that the excess number of atoms in the cluster should be 
a useful reaction coordinate. However, if one minimizes the model free energy with respect 
to R, w and po while holding AA^ constant, and then again solves perturbatively (now using 
e = (AA^)"^^'^ as the small parameter) the zeroth order central density is found to be given 
by 

du (poo) _ w (poo) - W (poo) ^^^^ 



^Poo Poo - P 
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oo 



which is physically incorrect since it implies that the final phase will not be the bulk liquid 
at the applied chemical potential. This happens because fixing the number of atoms in 
the cluster has the effect of changing the chemical potential as can easily be shown by 
formulating the problem with a Lagrange multiplier (see Appendix [U|) . 

A final possibility that suggests itself is to minimize at fixed equimolar radius. This 
avoids the issue of altering the chemical potential in the cluster. A perturbative solution, 
now using e = wo/Re as the small parameter, gives 

in 



An 



-RiAuj (poo) + 47ri?^7 



l + S 



Wo 

R 



O 



Wo 

R 



(42) 



with the zeroth and first order densities 
duj (poo) 



dpoo 
Poi 







(43) 



Uo 



Poo - poo (poo) 



^ Kjpoo) 
Ko 



Au (p 



'00 ; 



'00) 



2cJr 



Ko (P 



'00 j 



a width of 



Wo 



{poo - PooY 



-Ko (P: 



'00 ; 



2ujo - Aw (pooj 

and the coefficients for the expansion of the surface tension term are 



7 = Wo (2u;o - Aw (poo)) 

lujx- \Abj (poo) i^i (poo) 



(44) 



(45) 



Aw (p 



00 J 



+ 



K{p 



00) 



Poi 



2wo - Aw (poo) Ko (poo) 4 \2ujo - Aw (poo) Ko (poo) J Poo - Poo 

This energy-minimized path therefore gives the correct bulk density and seems the most 
likely candidate to be a good approximation to the nucleation pathway. 



C. Transition-state and steepest descent description 

The previous analysis shows that a description of the nucleation pathway in terms of 
minimum-energy configurations along some reaction coordinate is not trivial. Two candidate 
paths were found, one at fixed parameter R and another at fixed equimolar radius. This 
leads to the question as to which prescription is "correct" or whether there is a less arbitrary 
means of constructing such a description. 

The process of nucleation of a stable phase from an unstable one is conceptually similar 
to that of a chemical reaction or a structural transition in a finite cluster of molecules. 
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All of these involve the transition from a higher (free-) energy state to a lower energy 
state via an energy barrier. In principle, these processes should be described by dynamical 
theories. For the liquid-vapor transition, this would consist of a hydrodynamic description 
in which the free energy would enter via the local pressure. However, dynamical descriptions 
are computationally expensive and in the spirit of CNT, it is interesting to ask what can 
be learned simply from knowledge of the free energy functional governing the transition. 
For example, in CNT the free energy is a function of a single parameter, the radius, and 
the transition can be viewed as the growth of the radius from R = 0, the homogeneous 
metastable phase, to i? ^ oo, the homogeneous stable phase. The goal here is to generalize 
this picture for the case in which there are multiple parameters characterizing the transition, 
as is the case with the piecewise-linear density profile discussed above. It should be noted 
that the same question can be asked with regard to the SGA free energy functional Eq.Q, 
in which case there is a continuum of parameters, namely the values of p (r) for all points r. 
The standard approach to the descriptions between (meta-) stable states on an energy 



surface, widely used in the examples cited above[l3|, involves two parts. The first is the 
determination of transition states which are defined as saddle points in the free energy 
surface for which the Hessian of the (free-)energy function has a single negative eigenvalue. 
This is a generalization of the concept of the critical cluster in CNT. The second element is 
the determination of steepest descent pathways through the parameter space. These paths 
start at the transition state and consist of an initial small movement in the direction of the 
eigenvector with the negative eigenvalue. (Actually, there are two paths: one parallel to 
the eigenvector and the other anti-parallel.) The steepest descent paths are then the most 
efficient paths connecting the transition state to a local minimum - i.e. to a (meta-)stable 
state and, one expects, are closely related to the most likely paths when the processes is 
driven by thermal fluctuations. If a model for thermal fluctuations is available, then it is 
possible to define instead a most likely path 21| and of course, kinetic effects may alter the 
dynamics substantially. The steepest-descent paths simply represent the best guess of how 
the transition will proceed in the absence of dynamical information. 

In the simple case of liquid-vapor nucleation, one anticipates that there will be a single 
transition state (the critical cluster) and that the steepest descent path in one direction 
will lead to the homogeneous metastable phase while that in the other will lead to the 
homogeneous stable phase. In this case, it is the first part - the path connecting the initial 
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homogeneous metastable phase to the transition state - which is primarily of interest. 

For the sake of generahty, suppose that the free energy depends on n parameters and let 
r denote a particular set of those parameters so that for the piecewise linear model, n = 3 
and r = (po, R-i w). The transition state is a stationary point so must satisfy 

for all a = l,...,n. It can be efficiently located using so-called eigenvalue following 



techniques (22|. These are similar to gradient-following minimization techniques except 
that one minimizes in directions corresponding to positive eigenvalues of the Hessian and 
maximizes in directions corresponding to negative eigenvalues. 

The steepest descent paths are the quickest paths down the gradient starting at the 
transition states. The notion of "quickest" involves a measure of distance in parameter 
space which leads to the question of how to define a distance between two points Fi and F2 
given that the individual parameters will not in general even have the same units. Since the 
picture underlying the model is of a transition from one density profile to another, and since 
the concept of distance used to define the steepest descent directions should be independent 
of the model, it seems most natural to use the Euclidean distance in density space, 

d'[puP2] = J (pi(r)-p2(r)f rfr (47) 

which induces a distance measure in parameter space, 

d' (Fi, F2) = J{p (r; T,) - p (r; T,))' dr (48) 

It is clear that in the distance measure in parameter space is not Euclidean and in fact the 
induced metric is 

f dp (r; r) dp (r; T) ^ 
^7..(r) = y^^r ^dr. (49) 

The steepest descent path for a non-Euclidean geometry is then determined by 

dTb 1 dn 

ds Lab an an dVq 
y ava art, 



where s is the distance in parameter space ISj . As an illustration, the metric for the 



piecewise-linear spherical profile is given explicitly in Appendix |Dl The procedure is then to 
find the transition state which occurs at some point Fq, to make a small displacement of Fq 
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in the direction of the unstable eigenvector and then to use this as an initial condition for 
the solution of the steepest descent equations. Note that the metric, Eq. fHOl) . is only defined 
for density profiles which are at least continuous which is one reason that one cannot use 
the discontinuous zero-width profile implicitly assumed in CNT. 

Finally, it should be noted that there are alternatives to this procedure. Techniques such 



as the Nudged Elastic Band 23 



25| and the String Method 



are alternative, more heuristic. 



methods for determining steepest descent pathways. Both require a measure of distance in 
parameter space and so involve the same issues raised here. They are in general much less 
computationally demanding than the direct approach described above and have been used to 
determinepathways based on the inhomogeneous density, rather than a parameterization as 
used here jsl. I27I1 . The reason for using the present approach is that it seems most in keeping 
with the spirit of CNT plus the fact that with the simple piecewise-linear models, the direct 
integration of the steepest descent equations, Eq. flSUI) . is computationally straightforward. 



IV. COMPARISON OF THEORY AND SIMULATION 

To evaluate the SGA theory and the analytic approximations introduced above, a com- 
parison between the theory and simulations of planar and spherical liquid-vapor interfaces 
was carried out for a fluid with the Lennard- Jones pair potential 

„M=4.fn'^n'V (51) 



All forms of the theory require as input the bulk equation of state. In order to minimize errors 
arising from this external input, the empirical Lennard- Jones equation of state of Johnson, 



Zollweg and Gubbins (JZG)[28|] was used. Nevertheless, it is important to note that the 
JZG equation of state is based on simulations covering a limited range of temperatures and 
densities with corrections so as to describe the infinite-ranged Lennard- Jones potential while 
all of the simulation results used here are for a finite cutoff. Thus, in all cases, a mean-field 
correction was added to the JZG equation of state so as to account for the cutoff 



28 



It is 



therefore expected that the input is most reliable for large cutoffs and becomes increasingly 
unreliable as the cutoff is decreased. 

Figure [1] shows the free energy, density- averaged free energy and its first moment, 
Acu (p, p^) , Uq (p, py) and Cui (p, p„), for the Lennard- Jones potential with no cutoff for coex- 
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FIG. 1. The free energy and SGA coefficient and their zeroth and first-order density moments for 
a Lennard-Jones fluid with no cutoff at ksT = 0.75e = O.bSksTc and fi = ficoex- 

istence, fi = ficoex and at temperature ksT = 0.75e = O.SSfc^Tc which is just above the triple 
point (and where = 1.3e is the critical temperature). Also shown are K{p), Kq{p^p^) 
and Ki (p, p^). The free energy, Aw (p, p„), has two minima corresponding to the vapor and 
the liquid states separated by a barrier of about 0.36A;bT. The density-averaged free energy 
shows less variation and the first moment is about half as large. The SGA coefficient varies 
relatively little as a function of density thus showing that it is dominated by the density- 
independent potential contributions. As a consequence, the density-averaged value is nearly 
constant and the first moment is very nearly half as large, Ki (p, p„) ~ \Ko (p, Pv)-, over the 
whole range of densities. In FigJJl the same quantities are shown for = 1.2e = 0.92kBTc. 
In this case, the difference between the free energy moments is much less but the SGA co- 
efficient is still dominated by the density-independent contributions. 

A. Planar interface 

As a first test of the SGA, the profile and excess free energy of the planar liquid-vapor 
interface at coexistence was calculated for the Lennard-Jones potential truncated at various 
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FIG. 2. The same as FiglUbut for fc^T = 1.2e = 0.92/;:^^^. 



positions, r^, and shifted so that v{rc) = 0. Figure [3] shows the excess free energy per 
unit area, or surface tension, as a function of temperature from the calculations as well 
as that calculated in the piecewise linear model and values obtained from Monte Carlo 
simulations 29|. The accuracy of the calculations is evident. In fact, the accuracy of the SGA 
appears to rival that of the underlying DFT (see Ref.[l| for comparison). The piecewise- 
linear approximation always gives higher values of the surface tension, as it must since 
the SGA result is obtained via an unconstrained minimization of the free energy, but is 
nevertheless very close to the SGA result. 

Some calculated profiles are shown in FigJH and compared to simulation data reported 



in Ref. 29|. The SGA profiles are in reasonable agreement with the simulations although 
some discrepancy is apparent, particularly in the narrowest interfaces. The crudeness of the 
piecewise-linear approximation is apparent; even so, the widths of the interfaces are tracked 
reasonably well as a function of cutoff and temperature. 

To illustrate the convergence to the exact result as the number of links in the profile 
increases, calculations were performed while varying the number of links. The results are 
shown in FigjS] which shows that the simple single-link profile gives an error of about 5% 
and that this decreases to about 0.5% for 10 links. Linear extrapolation of the values as a 
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FIG. 3. The surface tension as a function of temperature for potentials cutoff at rc = 6cr (upper 
curves and symbols), rc = 4(T, Tc = Scr and Tc = 2.5a (lower curves and symbols). The symbols are 



the simulation data and error bars reported in Ref. 



29[, the full lines are from the SGA calculations 



and the dashed lines are from the piecewise-linear model. 



function of gives the exact value. 



B. Clusters 

Small clusters pose more of a challenge since there is no bulk region and most, if not all, 
of the atoms in the cluster are affected by the interface. Furthermore, all clusters are out 
of equilibrium except the critical cluster which is in a metastable state. The description of 
unstable clusters will be discussed below: here attention is focussed on the transition states: 
i.e., the critical clusters. 

Using the SGA, the properties of critical clusters were determined by solving Eq.(I7]) in a 



spherically symmetric geometry using a relaxation technique 



331 ]. Given an initial guess of 



the profile not too different from the critical cluster, this method relaxes to the critical cluster 



automatically. For the analytic model, the eigenvalue-fol 
was used. The protocol is the same as described in Ref. 

19 



owing technique described above 
1 

3[ including a small temperature 







10 20 

z* 



30 



FIG. 4. (Color online) Density profiles at the liquid-vapor interface calculated at different 
temperatures and values of the potential cutoff. From left to right, the curves correspond to 



ksT/e = T*= 0.7 and r* = 5.0, T* = 0.7 and r* = 2.5, T* = 0.8 and r* = 5.0, T* 



0.8 and 



the 
ref. 



2.5 and T* = 1.1 and r* = 5.0. The upper panel shows curves calculated in the SGA while 
ower one shows the piece-wise linear approximation. The symbols are the data reported in 



30l | and extracted from ref. 



31[ as the original is no longer available 



sal- 



correction to account for deficiencies in the equation of state for small cutoffs. 

Figure |6] shows the density profiles for critical clusters at two different values of the super- 



saturation as well as simulation data from ten Wolde and Frenkel 
Ref. 



34l ]. Note that following 



34j |. supersaturation is defined as the ratio of the vapor pressure to that at coexistence. 
The SGA is seen to give a good description of the critical clusters, although there are greater 
differences from simulation than in the case of the planar profiles. In particular, the SGA 
profiles have wider interfaces than occurs in the data and the larger profile, corresponding 
to lower supersaturation, has somewhat larger radius than indicated by simulation. Also 
shown are the analytic approximations which give a reasonable approximation to the SGA 
profiles but which are still wider and therefore compare less well to simulation. To put these 



results in context, the underlying MC-VDW DFT mode 
is in close agreement with the simulated cluster profiles 
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on which the present SGA is based 

I 
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FIG. 5. The relative error in the surface tension as calculated using piecewise-linear profiles with 
varying numbers of links, N, for r* = 6 and for T* = 0.7. The line shows the linear extrapolation 
of the profiles with 5,6,8 and 10 links demonstrating convergence to the exact value. 

Figure [7] shows the cluster size and excess free energy as a function of supersaturation as 
computed from the SGA, the piecewise-linear profiles, CNT and determined from simulation. 
The SGA gives a good estimate of the free energy barrier, but systematically underestimates 
the cluster size due to more rapid convergence to the bulk vapor density. The analytic profiles 
again give reasonable approximations to the SGA. As noted by ten Wolde and Frenkel 34 1 
CNT gives good estimates of the cluster size and poorer estimates for the barrier, especially 
at higher supersaturations. 



C. Nucleation Pathways 

The steepest descent pathways have been calculated by integrating Eq. flSOjl for several 
values of the supersaturation. Figure |H] shows the excess number of atoms in a cluster as a 
function of the distance along the pathway for two cases showing that AA^ is a monotonic 
function of the distance and therefore could serve as a reaction coordinate. Figure |9] illus- 
trates the actual pathways and shows some surprising features. Although the excess number 
of atoms is a monotonic function of the path, the size of the bulk region, parameterized by 
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FIG. 6. (Color online) Density profiles of the critical cluster for two different values of the su- 
persaturation for the Lennard-Jones potential with cutoff Tc = 2.5a and ksT/e = 0.741. The full 



symbols are the simulation data from Ref. [3j], the full lines were calculated using the SGA, the 
dashed lines are the piecewise-analytic approximation and the open circles are the result of the full 
MC-VDW0]. 

R, and the width of the interface, parameterized by w, are both non-monotonic functions 
along the steepest-descent pathway. The width in particular grows with growing droplet 
size for small droplets until it reaches a maximum for clusters of about 75 atoms and then 
slowly decreases thereafter. The size of the bulk region shows a local minimum at about 
the same point as the width has a maximum and is a decreasing function of cluster size for 
small clusters of 30 — 75 atoms. 

Figure M also shows the minimum energy pathways obtained by minimizing with respect 
to all parameters while holding AA^ fixed. The fact, noted above, that the inner density is 
incorrect for large AA^ is evident. However, the critical clusters are accurately determined 
since they are stationary points with respect to all parameters. This means that any energy 
minimized path, regardless of the constraints, will give the correct critical cluster. For 
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FIG. 7. (Color online) The properties of the critical cluster as a function of supersaturation for the 
Lennard-Jones potential with cutoff Tc = 2.5a and ksT/e = 0.741. The panel on the left shows 
the excess number of atoms in the cluster and that on the right shows the excess free energy. The 
squares are the simulation data from Ref. [sj], the full lines are the result of solving the SGA, 
Eq. (|19p . the dashed-lines are from the piecewise-linear approximation and, for comparison, the 
results obtained from the full MC-VDW DFTy] are shown as stars and the CNT result is shown 
as a dashed-dotted line. 



smaller clusters, the fixed-AA^ paths are in qualitative agreement with the steepest descent 
paths showing the same non-monotonic behavior of both the width and radius. However, 
rather than varying smoothly and continuously as the cluster size goes to zero, it was found 
that below a certain cluster size, the energy-minimized path jumped discontinuously to a 
solution consisting of a very low central density, only slightly larger than the gas, and a 
very large width. This is not an accident: for the small values of AA^, there are no energy- 
minimized clusters with liquid-like cores. Calculations of the energy-minimized path at fixed 
equimolar radius gave good agreement with the steepest-descent paths for large clusters, but 
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FIG. 8. (Color online) The excess number of atoms in a cluster as a function of distance along the 
steepest-descent paths where s = corresponds to the critical cluster. The figure shows the results 
for the Lennard-Jones potential with a cutoff at 2.5a for T = 0.683Tc for two different values of the 
supersaturation and illustrates the fact that the excess number varies monotonically with distance 
along the path so that it can sensibly be used as a reaction coordinate. 



show similar unphysical behavior (i.e. the absence of liquid-like cores) beginning at larger 
cluster sizes than for the fixed AA^ paths. The conclusion is that while energy-minimized 
paths can give qualitative behavior similar to the steepest descent paths, they do not give 



a reasonable physical description of the entire nuc 



eation pathway. This is somewhat at 



odds with the recent results of Ghosh and Ghosh [19|] who examine energy- minimized paths 
at fixed radius for sigmoidal profiles and who appear to obtain nontrivial minimizations at 
all values of the radius. Perhaps this is simply due the use of different equations of state 
and values of the squared-gradient coefficient or because the piece-wise linear model is too 
crude. Another possibility is that it is related to the fact that the sigmoidal profile does 
not enforce the boundary condition that dp{r)/dr = at r = as do the piecewise-linear 
profiles used here and so are less constrained. Note that this boundary condition is used 
when solving the SGA Euler-Lagrange equations, Eq. (fT9l) . for the critical cluster and that 
without it, derivatives such as dp{x, 0, 0)/dx do not exist at x = 0. 
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Another interesting feature of the steepest-descent pathways is that the the interior den- 
sity is that of the metastable vapor for very small clusters and increases rapidly as a function 
of cluster size until the cluster reaches about 100 atoms. Since the radius of the bulk region 
never fully goes to zero, this gives a very different picture of the formation of small clusters 
from that implied in CNT. Recall that in the CNT model, small clusters have small radii 
while the central density is always that of the bulk. Here, the picture is one of increasing 
density in a bulk region that is always of finite extent. Figure [TO] shows the equimolar radius 
which further emphasizes this difference. Whereas in the CNT model, the equimolar radius 
is equal to the radius of the cluster, and therefore goes to zero for small clusters, here it is a 
function of both the radius of the bulk region and the width and in fact is never small than 
about 2cr. This is very similar to results found using the NEB method and a much more 
sophisticated DFT modelj^, l^^]. Here, however, the physics behind this behavior is evident. 
As shown above, the surface free energy is proportional to the difference in densities inside 
and outside the droplet and inversely proportional to the width of the interfacial region. 
In order to minimize the free energy for small droplets, which is dominated by the surface 
tension contribution, as the size of the droplets decreases, the difference in densities must 
decrease as well. Since the free energy is inversely proportional to the width, the width 
stabilizes at a finite value and the AA^ — > limit is achieved via po ~^ Poo at finite w giving 
a non-zero equimolar radius. 



V. CONCLUSIONS 

The construction of the Squared-Gradient approximation starting with the Modified- 
Core van der Waals model for inhomogeneous fluids has been described. A relatively simple 
expression for the SGA coefficient was obtained which requires only the equation of state and 
interaction potential as input. The SGA model was shown to give quantitatively accurate 
surface free energies and density profiles for planar liquid- vapor interfaces of Lennard- Jones 
fluids as a function of temperature and potential cutoff. Similar comparisons were made 
for spherical clusters where it was found that the SGA was less accurate than the full DFT 
model nevertheless gives reasonable results. 

It was also shown that the SGA could further be approximated using piecewise-linear 
density profiles. This model is of course less accurate in the same examples (planar interfaces 
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FIG. 9. (Color online) The inner density, pQ, the radius, R, and the width, w, as a function of 
excess number of atoms in a cluster for the same conditions as described in Fig. [8l The full curves 
are the result of solving the steepest descent equations and the symbols mark the critical clusters. 
The dashed lines are the minimum-energy pathways calculated for fixed AA^. 



and spherical clusters) than the SGA but is not unreasonable given its simplicity. Aside from 
providing a rather simple means to explore the solution of the SGA, the main advantage 
of the piecewise-linear model is that it offers a simple bridge between DFT and the ideas 
behind Classical Nucleation Theory. 

Finally, the use of these tools to construct a description of homogeneous liquid-vapor nu- 
cleation was illustrated. The description was based on an analogy to chemical and structural 
transitions and made use of the transition-state/steepest-descent path methods used in those 
problems. Interesting non-classical behavior that was observed included the decrease of the 
interior density and the finite equimolar-radius as the size of the clusters tended to zero. 
Both effects were linked to the fact that the effective surface tension is density-dependent 
(unlike in CNT), a fact that is immediately evident in the case of the piecewise-linear ap- 
proximation but that would be harder to isolate in purely numerical calculations using the 
SGA or the original DFT. This density dependence is crucial in following small clusters as, 
otherwise, the surface tension would cause AQ to diverge as the interfacial width goes to 
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FIG. 10. (Color online) The equimolar radius (left panel) and the excess free energy (right panel) 
as a function of cluster size as calculated along the steepest descent paths. The calculations were 
performed for two values of the supersaturation and the critical cluster in each case is indicated 
with a circle. 

zero. This serves to illustrate the utility of such simple, yet not unrealistic, approximate 
methods in developing physical understanding of the more complex calculations. 
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Appendix A: Evaluation of the coefficients 

From Eq. f|T3|) . the explicit expression for the core-correction coefficients is 
&o = ^ [-^fHS (p) - (p) J + ScHs {d.- p, d) + ?>Pv id) + ^J^ Pv (r) r^dr (Al) 
hi = -f3v (d) - CHS {d^]P, d) - bo 
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The gradient theory coefficient is given by 

47r r , _ ^, 4^ 47r /"^ / rx 47r T 

K = —J CHs{r;p;d)rdr + — J (^'^o + Oi^ j - y 



/3f (r) r^rfr 



(A2) 



47r 
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c/f5 (j^; p; c^) ?"^c^?^ + ^ (660 + 56i 



47r 



I3v (r) r^cir 



— / C/f5 (r; p; d) r dr + — — (-/^/w {d) - chs [d-] p, d)) + -—bo - — (3v[r)r dr 
6 Jo 9 45 6 



(p) 



ap' 



■/(p) 



2 47r /""^ 

7rc/^/3t; id) + — / (2^^ - Sr^) (r) r^rfr 

45 30 ' 

Then, using 



dp' 



Chs {r; p; d) dr 



(A3) 



gives 



d' & 



poo 2 

A' = — / c«s (r; p; d) {6r^ - 2(1') rVr - -ird^c^s (rf-l P, <i) - ^j^fe. (p) (A4) 



(2rf2 _ 5r2) /3t; (r) r^dr 



For both the Percus-Yevick approximation and the more accurate White-Bear DFT, the 
hard-sphere DCF can be written as 



Chs {r; p; d) = ^^o + + as (^^) ^ 6 (d - r) 



(A5) 



givmg 



(r) 
dr 



dr 



Finally, 
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dp 
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where the isothermal compressibihty is 



= = -V 



d(5P 

'W 



So 



K 



180 IS'^ 45 ^ ^ 45 ^ ^ cir 



For packing fraction r] = ^pd^, the Percus-Yevick approximation gives 



as 



ri{l + 2riy 



2 

while in the White-Bear approximation, 

-3 + lOr/ - 15r/2 + 5r/3 31n(l-r/) 
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Appendix B: More general planar interfaces 



The piecewise-linear model is easily extended to include an arbitrary number of pieces. 
The density profile becomes 



P- 



z <0 



p{z)= < 



P~oo + (Pl - P-00) ^, < 2; < t^i 
+ {P 

P2 + (P3 - P2) 



Pl + (P2 - Pl) wi < 2 < u^i + u'2 



2— Wl— 11)2 
UI3 ' 



Poo, Wi + W2 + ... + W„ < 2; 



Defining 



this can be written as 



2;,; = > Wj 



P_oo, 2; < 



p(2) = < 



pi-1 + (pi - Pi-l) 



^i-1 < ^ < 



(Bl) 



Pcxd; ^ ^ 
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with the identifications po = P-ooj Pn = Poo- The free parameters are pi for 1 < i < n — 1 
and Wi for 1 < i < n giving a total of 2n — 1 parameters. Then, the excess energy is 



7 = = E jc^ (P (^)) - oo^oe. + (P (^)) (^^^) ' I (B2) 

i=l I * -^Pi-l 



which is simply the sum of the contribution of each link in the profile. 



Appendix C: Perturbative Solution for the spherical profile at constant particle 
number 

The free energy is 

- - f (-0 - (i) - - ^ (^0 - (i) - (i) 

(CI) 

which is to be minimized at constant particle number, 

iV = y (p (r) - Poo) dr = ^ (po - Poo) (2i? + w) {2Rw + 2R^ + w^) (C2) 
Introducing a Lagrange multiplier, A, and setting 



= — (Afi (i?, Po, w^) - A (iV - ^ (po - Poo) (2i? + u^) (2i?u; + 2R^ + 



2 I 2 
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for T = po,w and A gives, after some simplification, 



O.A. + 2^L..,g)+(-^(Z..X.(|))| ,C3) 



+A^(po-Poo) (3 + 3(1) 



r) 



3 dpo R \ dpo dpo \rJ dpo V-R/ 
^w {po- Poof / OKq ^ ^ dKi /w\ ^ 0X2 (w\'^ 



R 2w'^ \dpo dpo \RJ dpQ \R 



+Ai(,„-pj(6 + 8(|)+3^-j 
= A' - I (Po - fe) fi' (4 + 6 (I) + 4 ( 



fl) + (1) 



These are solved perturbatively using e = N as a small parameter where the expansion 
of the various quantities is assumed to take the form 

R = e-^Ro + Ri + eR2 + ... (C4) 

Po = Poo + epoi + •■■ 
w = Wo + ewi + ... 
A = Ao + eAi + ... 



Then, the lowest order equations are 



=UJ (poo) - (poo) + Ao (poo - Poo) (C5) 

du (poo) 

U — — 7: H Ao 

C'POO 

= UJ0 (poo, Poo) - Ko (Poo, Poo) + Ao J (poO - Poo) 

2wq 2 

= 1 - y (poo - Poo) Rl 
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giving 



Ao = _^(/^oo)-^(P-) (C6) 

Poo — Poo 

duj (poo) (poo) - UJ (poo) 



C^POO Poo — Poo 

^2 _ (Poo — Poo)^ 

° 2a;o(Poo,Poo) - (w(Poo) -'^(Poo)) ° 



° 47r (poo - Poo) 
Notice that the second equation can be written as 

df (poo) 



p - Ao (C7) 



5poo 

thus showing the shift of the chemical potential arising from the constraint. 



Appendix D: Metric 

The metric is calculated using Eq.f l49p . The result for a piecewise-linear profile with a 
single link is 



9pp = 


27r 
15 


(lOi?^ 




9pR = 


(Po 


- Poo) 


1 ^ (6i?2 + ARw + w"^) 




(Po 


-Poo) 


1 — (l0/?2 + lORw + 3^2) 
15 




(Po 




2 47r 3i?2 + + w"^ 


9rr = 


-Poo) 




1 




3 w 




(Po 




2Tt6R^ + 8Rw + 3w^ 


dRw = 


- Poo) 


1 




3 w 




(Po 




2 2tt lOR^ + 15Rw + Qw^ 


Qww 


-Poo) 




\ 

15 w 



[1] J. F. Lutsko, J. Chem. Phys. 128, 184711 (2008). 

[2] J. F. Lutsko, J. Laidet, and P. Grosfils, J. Phys.: Cond. Matt. 22, 035101 (2010). 

[3] J. F. Lutsko, J. Chem. Phys. 129, 244501 (2008). 

[4] J. S. Rowlinson, J. Stat. Phys. 20, 197 (1979). 

32 



[5] J. D. van der Waals, Z. Phys. Chem. 13, 657 (1894). 

[6] V. Ginzburg and L. Landau, Zh. Eksp. Teor. Fiz. 20, 1064 (1950). 

[7] J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 (1958). 

[8] R. Evans, Adv. Phys. 28, 143 (1979). 

[9] H. Lowen, T. Beier, and H. Wagner, Europhys. Lett. 9, 791 (1989). 

[10] H. Lowen, T. Beier, and H. Wagner, Z. Phys. B 79, 109 (1990). 

[11] J. F. Lutsko, Physica A 366, 229 (2006). 

[12] P. M. W. Cornelisse, C. J. Peters, and J. de Swaan Arons, J. Chem. Phys. 106, 9820 (1997). 

[13] D. Wales, Energy Landscapes (Cambridge University Press, Cambridge, 2003). 

[14] J. -P. Hansen and L McDonald, Theory of Simple Liquids (Academic Press, San Diego, Ca, 
1986). 

[15] R. Evans, in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson (Marcel Dekker, 

New York, New York, 1992). 

[16] J. F. Lutsko, J. Chem. Phys. 127, 054701 (2007). 

[17] J. A. Barker and D. Henderson, J. Chem. Phys. 47, 4714 (1967). 

[18] M. M. T. da Gama and R. Evans, Mol. Phys. 38, 367 (1979). 

[19] S. Ghosh and S. K. Ghosh, J. Chem. Phys. 134, 024502 (2011). 

[20] D. Kashchiev, Nucleation : basic theory with applications (Butterworth-Heinemann, Oxford, 
2000). 

[21] M. Heymann and E. V. Eijnden, Phys. Rev. Lett. 100 (2008). 

[22] J. C. Mauro, R. J. Loucks, and J. Balakrishnan, J. Phys. Chem. A 109, 9578 (2005). 

[23] H. Jonsson, G. Mills, and G. K. Schenter, Surface Science 324, 305 (February 1995). 

[24] G. Henkelman and H. Jonsson, J. Chem. Phys. 113, 9978 (2000). 

[25] G. Henkelman, B. P. Uberuaga, and H. Jonsson, J. Chem. Phys. 113, 9901 (2000). 

[26] E. Weinan, W. Ren, and E. V. Eijnden, J. Chem. Phys. 126 (2007). 

[27] J. F. Lutsko, Europhys. Lett. 83, 46007 (2008). 

[28] J. K. Johnson, J. A. ZoUweg, and K. E. Gubbins, Mol. Phys. 78, 591 (1993). 

[29] P. Grosfils and J. F. Lutsko, J. Chem. Phys. 130, 054703 (2009). 

[30] M. Mecke, J. Winkelmann, and J. Fischer, J. Chem. Phys. 107, 9264 (1997). 

[31] K. Katsov and J. D. Weeks, J. Phys. Chem. B 106, 8429 (2002). 

[32] J. Fischer, private communication (2007). 

33 



[33] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C 

(Claredon Press, Oxford, 1993). 
[34] P. Rien ten Wolde and D. Frenkel, J. Chem. Phys. 109, 9901 (1998). 



34 



